

    ### Starting value

    KStart = function(eX,Z,N){
        
        mie = rank(eX)
        
        Z.order         = Z[order(eX)]
        Z.gaps          = diff(Z.order)
        Z.change        = c(0,which(Z.gaps!=0),N)
        Z.const.length  = diff(Z.change)
        Z.const.argmax  = which.max(Z.const.length)
        Z.const.maxlength = max(Z.const.length)
        
        Z.const.endpoints = Z.change[c(Z.const.argmax,Z.const.argmax+1)]
        
        K.start = round(2*N/Z.const.maxlength)
        
        return(list(K.start=K.start, Z.const.endpoints = Z.const.endpoints))
            
    }

    
    ### Maximal number of subclasses

    KMax = function(eX,Z,K.start,Z.const.endpoints=c(0,0)){
    
        K = K.start
        upper = Z.const.endpoints[2];    lower = Z.const.endpoints[1]
        
        flag = 0
        while(flag == 0){
            
            if( abs(ceiling((upper+1)*K/N)-ceiling(lower*K/N)) < 1.5){
              q    = unique(quantile(eX,prob=(0:K)/K))
              q[1] = q[1]/2; q[length(q)] = 1
              K = length(q)-1
              ki   = as.numeric(cut(eX,q))
                if(min(table(Z,ki))>0) flag = 1
            }
            
            K = K - 1
        }
        
        return(K+1)
        
    }
    
    
    ### Number of subclasses based on subclassification method 
    
    KEst = function(Z,eX){
        
        Z = Z[!is.na(eX)];  eX = eX[!is.na(eX)]
        N       = length(Z)
        k.start = KStart(eX,Z,N)   # Find starting value for K
        k.max   = KMax(eX,Z,k.start$K.start,k.start$Z.const.endpoints)
        
        return(k.max)
        
    }
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    